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Abstract 

The mean field approximation is numerically validated in the bosonic case by considering the time 
evolution of quantum states and their associated reduced density matrices by many-body Schrodinger 
dynamics. The model phase-space is finite-dimensional. The results are illustrated with numerical 
simulations of the evolution of quantum states according to the time, the number of the particles, and 
the dimension of the phase-space. 
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1 Introduction 

The mean field approximation is known to be a good way to approximate the many-body Schrodinger 
dynamics when the number of particles is large enough (see [2, 6, 9, 10, 15, 19, 20, 21, 22, 24, 25, 26, 32, 37, 
43, 46, 29, 30, 33, 50]). 

It consists in looking for the solutions to the non-linear Schrodinger equation for one particle called the 
Hartree equation. We are interested in the density matrix associated with the wave function, this matrix 
satisfies the quantum Liouville equation dual to the Von Neumann equation.The partial trace operators of 
this matrix, called the reduced density matrices, satisfy a hierarchy of equations. For instance, by considering 
the case where the initial state for the Schrodinger equation is a Hartree ansatz(a product state) which is 
suitable for a bosons condensate, the limit, when the number of particles N goes to the infinity of these 
matrices converge in trace norm to the product of the density matrix associated with the solution to the 
Hartree equation. And this asymptotic density matrix satisfies the time dependent Hartree equation [8]. 
When the particles are bosons, the suitable space for the bosons is the symmetric Fock space on the phase- 
space. Moreover for the sake of numerical computations, a finite-dimensional phase-space will be used instead 
of an usual phase-space of type L 2 (R d ). So here the phase-space will be Z = £ 2 ({0, • • • , K }) ~ C K where K 
is a given integer representing the number of sites. Each particle can live in one of the K sites. 

For the numerical implementation, an explicit basis of the iV-fold sector of the Fock bosonic spaces is 
specified. This basis allows the numerical computation of the full iV-body quantum problem for N large 
enough to validate various mean field regimes, in spite of a rapidly increasing complexity. 

The resolution of the 7V-particles Schrodinger equation will rely on a splitting method, one part for the free 
Hamiltonian and the other one for the two particles interaction term. 

For the simulations, the considered real bounded potential associated with the interaction term will be V 
defined on Z/ATZ by V(i) = j|j- if * ^ 0 and V(0) = 0. 

According to previous results related to the propagation of the Wigner measures [3, 4, 5, 6] knowing the 
Wigner measure at time t = 0 determines the Wigner measure at time t and all asymptotic reduced density 
matrices. For many examples, like Hermite states, twin Fock states or states studied in quantum information 
theory (see [1]) their Wigner measure as well as the order of convergence of reduced density matrices is known 
explicitly. The evolution of the Wigner measure, and consequently of the asymptotic density matrices, is 
evaluated after integrating numerically the mean field non linear Hartree time-dependent equation. In order 

*IRMAR, Universite de Rennes I, UMR-CNRS 6625, Campus de Beaulieu, 35042 Rennes Cedex, France 
'Fak. Mathematik, Univ. Wien, Oskar-Morgenstern-Platz 1, 1090 Wien, Austria 


1 



to preserve numerically quadratic quantities like the symplectic form, the latter is solved with a symplectic 
4 th order Runge-Kutta method ([35]). 

To estimate numerically the error of convergence of the reduced density matrices in the mean field limit, 
a discretization of a time interval [ 0 , t max ] is considered, in the examples t max = 1 is chosen, then the 


quantity max f£ [ 0itraiii ] 7/^ (t) ~ 7 oo J (t) is observed. Here 7 ^’(t) denotes the time-evolved p particles 


„(p)/ 


(p)/ 


reduced density matrix for N bosons, while 7 ^ (t) is its theoretical limit when N goes to the infinity, and 
||.|| 1 denoting the trace norm. 

For the evaluation of the order of convergence, the logarithm of the previous error estimate is drawn according 
to log(N). This gives a straight line whose the slope is the order of convergence in 1/N. 

These numerical results agree very well and illustrate the theoretical analysis carried out in [1], 

By increasing I \, we wish to approach a continuous model. The complexity of the computations increase in 
the same time that K and N increase because of the dimension of the 7V-particles bosons Fock space on C K 
which is a binomial coefficient. 


2 Framework 

The bosons Fock space on an Hilbert space Z is defined as V S {Z) = ©„> 0 V" % where \J n Z is the symmetric 
n-fold Hilbertian tensor product of Z which is the range of the projection defined on the Hilbert tensor 
product Z® n , by: 


S n {£. 1 ® £2 ® ... <g> £n) = — ^ £,7(1) <8> £<t( 2 ) ® ... ® Ccr(n), 
where is in Z for each i in [1, n] and E„ is the set of the permutations of n elements. 

For z in Z and e positive, the e-scaled annihilation and creation operators are defined for all <f> in Z and 
n in N by: 


a(z)4>® n = V^(z 

a*(z)4>® n = ^/e(n + l)S n+1 (\z) ® $® n ). 
These operators are then extended by linearity and density to \J" Z. 


These operators satisfy the canonical commutation relations (CCR): 


[a(zi),a*(z 2 )\ = e{zi,z 2 )Id , 

[a(zi),a(z 2 )] = 0 , 

[a*(zi), a*(z 2 )\ = 0 . 


( 1 ) 

( 2 ) 

(3) 


The second quantization of an operator A £ C{Z) or a self-adjoint operator (A,D(A)) in Z is defined 
by: 

n 

dT{A)^, a i aD(A) =eY J Id® l ~ 1 ® A® Id® n ~ l . 

*=1 

The second quantization of Idz is the number operator: 


N| V "2 = enldynz . 


2.1 Orthogonal basis of the iV-fold sector 

Use the following notations: Z x = Z/R'Z . 

For a = («!,••• ,cvk) in N^, the length of a is written |a| = aq + • • • + an and the factorial of a 
a! = a:i! • • • a'! . 

Let (ei, • • • , ex) be an orthonormal basis of C K . 

Set Z = <C K . Then an orthonormal basis of \/ N Z can be built from this basis which is labelled by the 
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multi-indices a in such that |a| = N. 

With the creation operators, an orthonormal basis can be written as: 


e„ := tCSllID := 


1 


V d Q la! V £ , I“Iq'! 

where |Q) = (1,0, 0, 0,...) is the vacuum of the Fock space. 


a*(e 1 ) ai ---a*(e^) a *|f2), 


Then the dimension of \J N Z is card({a £ /|a:| = N}) = 


N + K-l 
I< - 1 


And card({a £ N^/la] < N}) = ( N ^ K ) is the dimension of ®^ =0 V" Z ■ 


2.2 Hamiltonian 

As a binomial number, the dimension of the N particles bosonic sector increases rapidly but not too much, 
as N increases (see Table 6.4 for numerical values).The complexity has to be handled carefully if we want 
to approach the mean field limit by taking N large or the continuous model by taking K large. 

Define A k the discrete Laplacian operator on C K by: 

Vz £ C K , Vi £ Z/ifZ, (A K z)i = z i+ \ + Zi-i . 

And let Hq = dT(—Afc) be the free Hamiltonian. 

The interaction term denoted by V equals: 

V=^ Vijd*(ei)a*(ej)a(ei)a(ej) , 

(i,j)£Z 2 K 

where E? = Vp = V{i — j). 

In this framework changing the value of V (0) add an irrelevant phase factor in the time evolved wave function. 
In the sequel V(0) =0 is assumed. 

Or as a Wick quantized operator (14): 

v = ® e i)(o ® e i\) z ® 2 ) Wick . 

(ij)e z 2 k 

The considered linear Schrodinger equation is: 

i ed t V = H t A> , (4) 

where the complete Hamiltonian is defined on the bosonic Fock space by : 

H e = dT(- A k ) + V . 

3 Finite dimensional mean field equation 

3.1 Energy of the Hamiltonian 

The energy of the Hamiltonian corresponds to the symbol of the complete Hamiltonian: 

H{z,z) = (z,-A k z) + lj2 V ij\zi\ 2 \zj\ 2 , 

while recalling our convention Vu = E(0) = 0 for all i in 'Lk- 
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3.2 Hartree equation 

The mean field equation in Z is written as: 


i d t z = d s H{z , z) . 

For each component in C K we obtain: 

idtZi = d Si H = W ~ - 2 N 2 + \ Vi'j\zi'\ 2 \zj\ 2 ) 

V i'^j 

= dzi(52{zi' - zf-i){zi / - z v - 1 ) - 2 ZiZi + Vi-jZi'Zi-ZjZj) 
= - Zi- 1 - (^+1 - Zj) - 2z, + ^ V^lz,] 2 

= — (^i+l + ^i-l) + 

ii = -i[-(A*rz)i + (Xl^l^l 2 )^] 

By writing 2 = g + ip where q and p belong to it becomes: 

Qi = -(p i+ i +Pi-l) + ^2 +Pj)Pi 

Pi = -[-(Qi+i+Qi-i) + ^2Vij{Qj+P 2 j)Qi) 

j¥=i 


A 4 th or 6 th order Gauss RK method is used by using the coefficients given in [35] to solve the Hartree 
equation. 

A symplectic method is used to preserve the quadratic part of the energy, the symplectic form and the phase 
space volume. 

3.3 Wigner measures 

For / G Z the field operator is defined by 4>(/) = (a*(/) + a(/)) which is essentially self-adjoint on 

r MZ) = @nixV n z 

The Weyl operator is defined by W(f) = e l9 ^\ 

Let {g E ) Ee e be a family of normal states on F z (Z) with £ C (0, +oo), Oef. 

A measure p is a Wigner measure for this family, p € A i(g e ,£ G £), if there exists £' C £, 0 G £' such that 


v/gz 


lim Tr 

,£—>-0 


g E W{V2nf) 


0 2i7rRe {f,z) 


dp(z) , 


see [42]. 

The following result valid for separable Hilbert spaces Z, apply to our finite dimensional Z ~ C K . 

Theorem 3.1 [3]. If (g E ) E ^£ satisfies the uniform estimate Tr^N 5 ] < Cg < +oo for some S > 0 fixed, 
A 4(g e , £ G £) is not empty and made of Borel probability measures (Z separable) such that J z |z| 2<5 G?/i(z) < Cg. 

For each p in N, the reduced density matrix associated with a state g E is a trace class operator in C 1 (\J P Z) 
defined by the duality relation: 


Tr 



Tr [g e b Wick ] 

Tr [g e (\z\ 2 P) Wick ] 


(5) 
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where be C(\J P Z). 

The asymptotic reduced density matrix associated with the Wigner measure g equals: 

/ z \z®<-)(z®r\d»(z) 

7 “ = JzWM’) ' (6 

In hnite dimension, if the family (g e ) e ^£ satisfies M(g e ,e £ £) = {g} then the (PJ)-condition (see [3]): 


Vp£ N, lim Tr[p e N p ] = / |z| 2p dg(z) 
ee£,e —>0 J z 

is always satisfied. 


3.4 Reduced density matrices 

Theorem 3.2 [5]. If the family (g e )e&£ satisfies M(g e ,e £ £) = {p} with the (PI)-condition: 


Vp £ N, lim Tr[p e N p ]=/ |z| 2p dg(z) , 

J z 

then Tr [g e b Wlck \ converges to J z b(z) dg(z) for all polynomial b(z) and 


J im Jhf-Toll C 1 =0 


for all p £ N. 

Theorem 3.3 [5, 6, f2]. Assume M(g e ,s £ (0, e) = {po} and the condition (PI) ■ Then M(e~' l ^ He g s e z ^ Hs ,e £ 
(0,e)) = {/-<*}• The measure g t = , t > (f,0)*po is the push-forward measure of the initial measure go where 
<E>(i, 0) is the hamiltonian flow associated with the equation 

id t z k (t ) = -A K z k (t) + ^2 v kj\zj\ 2 z k . (7) 

j 

After propagation of the Wigner measures, for any p £ N, the convergence of the reduced densiy matrices 
is obtained at any time t: 

J z \z® p )(z®*>\dg t (z) 

%U f z \z\*PdM*) 

Theorem 3.4 [1]. Let (a(n)) n 6 N* be a sequence of positive numbers with lima(n) = oo and such that 
(~|y)neN* is bounded. Let (g n )ne N* and ('Yx)pe n* be two sequences of density matrices with g n £ _Sf’ 1 (\/" Z) 
and 7^1 £ JZ >x (f\J v Z) for each n,p £ N*. Assume that there exist Co > 0, C > 2 and 7 > 1 such that for all 
n,p £ N* with n > 7 p: 



ryiP) _ ry(p) 
In loo 


<C 0 - 


C’P 


1 a(n) 


Then for any T > 0 there exists Ct > 0 such that for all t £ [—T, T] and all n,p £ N* with n > 7 p, 

C p 


T^\t)-T^(t) 


< Ct- , x , 
1 a(n) 


( 8 ) 

(9) 


where 


7 &>(*)= Jjz® p )(z®*>\dg t (z), 

with pt = ($t)iPo is the push-forward of the initial measure go by the well defined and continuous Hartree 
flow <E> t on Z. 
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4 Numerical methods 


4.1 Method to solve the Hartree equation 

To solve the mean field equation (7), a Runge-Kutta method is used. 

Let bi, dij (i,j = 1, • • •, s ) be real numbers and Ci = ]Cj=i a ij- 
An s-stage Runge-Kutta method with a time step h to solve a first-order ordinary equation y' = f(t 7 y) , 
2/(to) = Ho is given by: 

S 

h = /(t 0 +Cjh,y 0 +hy^a ij k j ) ,i = l,...,s 

3 =1 

S 

2/i = 2/o + hy^bjkj 

i —1 

represented as: 


Cl 

an 

au 

Cs 

a s i 

• Clss 


bi ■ 

■ b s 


Here the system is autonomous, and according to [35] the coefficients used for the Gauss RK method are: 


0 

1/2 

1/2 

1/2 

0 

1/2 

0 

1/3 

2/3 

1/3 

-1/3 

1 

1/2-V3/6 
or 1/2 + -\/3/6 

1/4 

1/4 + v^/6 

1/4 — a/3/6 
1/4 

1 

0 

0 1 

1 

1 

-1 1 


1/2 

1/2 


1/6 

2/6 2/6 1/6 


1/8 

3/8 3/8 

1/8 




In our case, the function / corresponds to /(z) = —i A K%k + Ylj Vkj\zj\ 2 Zk 

As a function in by replacing z by q + ip, 


f(<LP) = fqiQiP) + i f P {q,p) 
fqMiP ) = ~{Pi +1 +Pi- 1) + Vi i( q J + P 2 q)Pi 

f Pi (q,P ) = Qi +1 + Qi -1 - ^2 Vi i( q j +P 2 j)<li 

For an implicit Runge-Kutta method, a Newton method is applied to find the coefficients ki for each step 

of the RK method to the function g yo : (fci)j=i ... g (A;* — f(yo + h a ijkj) ) 

V ' J 

Given the time step h small enough, the starting point of the Newton method is chosen by setting 
h = f(yo) for all i. 


To apply the Newton’s method the differential of / is computed by using the following partial derivatives 
of/: 


dfq t 

dqk 

dfpi 

dqk 

df q , 

dpk 


— (1 ^i 7 k)^^ikqkPi 

= &i+l,k T 3i— l,k T (&i,k 1 )2Vikqkqi Si,k ^ ( K;y ( y Qj A Pj ) 

= ~{Si+l,k + Si-1,k) + (1 — Si,k)2VikPkPi + Si,k ^2 ^ij +?*|) 


df Pl 

dpk 


(1 Si^k)^^ikPkqi 
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Then the differential of g Va is: 


^9yo ((^'*)i=i,-" ,s) — I hauDf(yo T h ^ kj ) 

i=i 


i,Z=l, 


where g Vo is considered as a function from Kr As . 


4.2 Resolution of the Schrodinger equation in \J A Z 

4.2.1 Composition method 

For a given U/ in \j N Z, the full TV-body evolved state is computed in the basis (e a )| a | = jv. After 

writing d' = ^ a s-a a modified splitting method for which the numerical error is carefully controlled 

(see 5), involves only multiplications by the diagonal matrix e~ lj p v and the sparse matrix <£T(—A^). 

In order to handle the high complexity of the problem (see table 6.4) no matrix, but only vectors or the 
sparse matrices dT(—Ajf) and the matrix (V) ? ),;j are stored. 

The complete evolution e~ l i He is computed by a composition method based on the Strang splitting method: 


e l e He = lim (e * 2E » v e l ep H °e l ^p v V 

p—too v 7 


The 4 t/l order composition method is given by: 


p—too ' ' 


where the coefficients of the method are satisfying the two equations (see [35]): 


and are given by: 


# 1+82 + ®3 — 1 

Clf + U2 4“ ^3 = 0 


1 2 1 / 3 
ai ~ ° 3 ~ 2 - 2 1 / 3 ’ “ 2 “ 2 - 2 x /3 ■ 


( 10 ) 

( 11 ) 


( 12 ) 


4.2.2 Computation of the free evolution e - 1 j dr (~ A K) 

The numerical computation of e~' l i dr (~ A ^) = r(e itAlsr ), relies on the following two remarks: 

• the dimension of the TV-fold sectors prevents the storage of any square matrix. 

• the matrix of T(e ztAjf ) is actually non trivial sparse matrix in the basis (e a ). 

The matrix of Ax is given by: 

1 0 ••• 0 1 \ 

•. 0 

•. ■•. ■•. 0 

1 

0 ••• 0 10/ 
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A k = 


/ 0 
1 

0 


0 
V 1 






We are interested in the matrix of the second quantization of the discrete Laplacian on the basis of the 
bosons space to implement it numerically as a sparse matrix containing only 2K( N ^ < K f 2 'j elements whereas 

a full matrix contains . 

Then e -*^ dr (- A x) w ;p be computed at each time step by a 4f h order Taylor expansion. 

This expansion is then replaced in the composition method. 

For an operator A : C K — > C K , A = (Aij)ij , 

K 

dT(A) | V "c K = X! A hi a *( e i) a ( e j) ■ 

*>i=i 

This yields: 

K 

dT(-A k ) = -J2 a *( e j+i) a (ej) + a *( e j)a(e j+1 ) . (13) 

j =i 

Lemma 4.1 For all multi-indices 7 and a the following equality holds: 

a(e) 7 a*(e)“|Q) = • 

(a- 7 )! 

Proof. 

a(e) 7 a*(e)“ = a(ei ) 71 ... a(eif) 7 x a*(ei) Ql ... a*(eif)“ K 

K 

= a(ei) 7 i a*(ej) ai which is a commutative product because of CCR (1) 
i=l 
K 

= (£) a(ei) Ji a*(ei) ai 

i—1 

by using the following separation of the variables: T(Z) = T(Cei) ® ... <g) T(Ceif) ■ 

In this space let |Q) be |fii) ® ... < 8 > |n^). 

Let us consider 7 * > 1 and ai > 1, 

a(e i )' 1i a*(e i )°‘ i \n i ) = a(e i ) 7 i_ 1 a(e i )a*(e i ) ai |fl l ) 

= a(e i ) 7 i- 1 a*(e i )“ i a(e i )|fl l ) + a(e i ) li ~ 1 {a(e i ), a*(ej) ai ]|ftj) 

= £a l a(e i ) 7 i- 1 a*(e i )“ i ~ 1 |f 2 i) . 


By induction, we obtain 

a{ei) li a*{ei) ai \^li) = one x (on — l)e x ... x 2e x ea(ei) Ji ~ ai |fij) = 0 , when a* < 7 \ 

and 

a(e i ) 7 i a*(e i )“ i |fl i ) = a;£ x (a* - l)e x ... x (a* - ( 7 * - l))£a*(ej)“ i_7i |fi*) 

= e 7i 7—— —rra*(ej) a<-7i |f 2 j) ,when 7 * < a*. 

(a* - 7*)! 
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The above separation of variables leads under the condition 7 < a to: 
a(e) 7 a*(e)“|fl) = (gi^ 1 a(e i ) 7 i a*(e i ) ai |fb) 

e 7i a,;! 

(a* -7i) ! ' 


£ ' 7 ' (^)T (n (|fii)®...«|njf)) 


= el 7 l - 


(a - 7 )! 


^(e)“- 7 |fi) . 


Proposition 4.2 Tor all multi-indices a and /?, f/ie matrix elements of dT(—A k) are given by: 

dr(- A K )a,p = -sYX5+_ e ^ a _ e _ \/Pi{(di +1 + 1) + fip-ei+i,<x-eiV Pi+l(Pi + 1)) > 

i 

where 5+ p = (5 Q , ;/ gl N i<r (a) for a and (3 multi-indices in Z K . 

Proof. 

According to (13) and to Lemma 4.1, we obtain: 


□ 


a*(e i+ i)a(e i )a*(e) a |f 2 ) =S ei < a £ 


a: 


,a*(e) 


ei-l-l+Q! — €i 


(a - e-i)! 
=S 1 < ai £a i a*(e) ei+1+a ~ ei |fi) 


|fi> 


and 


l ‘(e i )a(e i+ i)a*(e)“|fl) =S e < a e 


v! 


K (e) 


ei~\-OL— ei-|-i 


ei+1 - (c«_e i+1 )!' 

=di<a i+1 £oii+ia*(e) ei+a ~ ei+1 |fl) . 


|fl) 


Then 


(e a ,dT(- A K )ep) = - (^e a , ^E a*(e i+ i)a(e,;) + a*(ei)a(e i+ i)^ e^j 


— £ 


- E ( a *( e ) Q ^l {h<^*{^ +ei+x - ei + S 1 <p i+ J i+1 a*(e)^- e ^) Si) 

= £ n J^W. E^-^-e.+^xMl - e t + e i+ i)\ 

+ Sp_ ei+i , a _ ei /3i + i£ N ^/a\(/3 - ej+i + e*)!) 

= ,,/W(ff ~ e » + e ^+l) ! + &+1 \/(/3 ~ e i+ i + ej)!) 


= EE s »- 


— ei+i,o:—e 






/ A +1 

/3i+l 


dr(-Aic) Q>(S — ~ g E(^-e. .a-e^ 1 A(&+1 + 1) + .a-e, V^i+l(A + !)) ■ 
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□ 

Numerically only the indices of the multi-indices a and [3 corresponding to the nonzero components of 
dr(— Ak) with their values, are stored in an array. 

In the algorithm instead of running over the multi-indices a or (3 with a length TV, the multi-indices f}' with 
a length N — 1 are run over. And for each % in [1, A'], the changes of multi-indices /3' = (3 — e t or (3' = (3 — e i+1 
are used, then the indices of the corresponding multi-indices a and /3 with length N are looked for. 
Therefore an array composed of 2 K ( W J A ^ 2 ) triplets of elements is numerically stored. 


4.2.3 Computation of the interaction factor e ’* v 

Denote a# = aft{e.j) and IV, = a*cti . 

By using the relations CCR (1): 

a*a*aittj = a*(dia* — eSij)aj = a*aia*aj — edija*aj 
= NiNj - eS^N, = N^Nj - e6 XJ ) . 

Then V can be rewritten as: 

V = - y] VijNi(Nj — e5ij) . 

(i,3)£%K 


And since Aje a = eaie a then V is diagonal in the basis (e Q ) Q : 


126 a — I ^ A \ 1 '{ j - O / (- Oy -Ay ) 

{ij)€% 2 K 


— I 9 Vijai(ctj 5ij) 




and 


-i-E-V 
e £ e n = e 




= Vijouaj+Y^ i&K Vuai(ai- 1)) 


4.3 Numerical computation of the reduced density matrices 

Consider b Wtck = a*(e) s a(e) 7 with |d| = lyj and its associated homogeneous polynomial: 

b(z) = z 6 z 7 = 4 1 ■ ■ ■ z S kzT ■ ■ ■ Zk • 

Let us compute the quantity Tr (g E b Wtck ) when g e is a normal state. Using an orthonormal basis of the 
TV-fold sector \/ A Z, g e is a linear combination of operators |$)(4 , |. It suffices to compute Tr(|$) (ty\b Wlck ) = 


Lemma 4.3 Set b w%ck = a*(e) 5 a(e) 7 with |<5| = |y| and let $ and 4^ be in \J Z then: 


(^,b Wick ^) =s H ^ 4v +(5 $ a , +7 

|a'|=]V-|5| 


\J ( a' + <5)!(a' + 7 )! 


in the orthonormal basis ( a , ^ of\J N Z. 

VV?^! 1 'J\ a \=N 
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Proof. Given ’I' and d> in the N-particles bosons space and the formula 4.1 


a*(e)' 5 a(e) 7 a*(e)“|f 2 ) = 8 1 < a e 


M. 


(a - 7 )! 


1 *(e)' 5 +a_ 7 |fl) 


we can write : 


= ( E • £h1 E *« ? w7fr7T“‘( e > i+ '’"'l n > 

\ot\=N Vfc “• |/3|—AT,/3>7 U ' 


-ItI 


VW- 


e N Z Z 

|q|=AT |/3|=iV,/3>7 r;. 

v / (o + 7-(5)!^| 
a+7—<5 /—-. 


a*(e)“fl, a*(e) < 5 +/ 3 _ 7 |fl) 


— c-It'I 


= £"' Y 


|a|=AT,Q!>(5 


= £ W Y ^a$a+ 7 -5 


|a|=A^,o:><5 


v a!(a — <5)! 

A/cd(a + 7 — <5)! 
(a — <5)! 


= £ l71 Y ^a'+5^a'+7 

|a'|=JV-|5| 


y/ («' + <5)!(a' + 7 )! 


because 


(a*(e)“fl, a*(e) <s+/ 3 - 7 |f 2 ) ^ 0 


if and only if a = 8 + /3 — 7 SO /3 = a + 7 — <5 and ft > 7 means a — (5 > 0 . 

The last line is obtained by a change of multi-indices by setting for each a, a' = a — 8 because a > 8 , and 
then \a’\ = |a| — |<5| = N — |<5|. □ 

Numerically, all multi-indices of N A with length not larger than a given N rnax are stored in the lexico¬ 
graphic order. 

For our algorithms, we pay attention to preserve this lexicographic order (or reverse). 

For a given N < N max , the list of relevant multi-indices (with length N) is extracted and handled in the 
lexicographic order. 

For a given 8 , numerically the above summation is performed over multi-indices cd such that \a'\ = N — 1<5| 
in the lexicographic order. 

Then for each a', the multi-indices a of length N written as a = a! + 8 are looked for. These a are exactly 
the multi-indices such that a > 8 and |a| = N. 

Note in particular that the mapping a' 1 —> a' + 8 preserves the lexicographic order. 

First let us compute the matrix elements of yj? in the orthonormal basis (e a ) a . 

The matrix element corresponding to the line (3 and column a is: 


(<n 




a*(e)“ 


le)* 




VePcd /\vWI 
Tr ( 0 e (t)b wick ) 
£PN(N - 1)... {N - p + 1) 


according to the duality relation (5) of the reduced density matrices with b = 
b(z) = {z®p, bz®P) € V p , p . 


_ 

tfa! “/\ \A p P ! 


a *( e )° 


n 


If 6 = 




e p a\ / \ %/e p /3! 


its Wick quantized is : 


b Wlck = —/==a*(e) a a(e)P . 
V«!p! 


n 


and 
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And then 


7?(*)(&«) = 


p! Tr(p £ (t)a*(e)“a(e) /3 ) 


Vap £PiV(IV - 1)... {N - p + 1) ' 

Then owing to Lemma 4.3, all the elements of the matrices can be numerically computed. 

In the case where the initial state is a Hermite state g e = |z 0 JV )( 2 :® JV |, z® N needs to be expanded in the 
orthonormal basis (e a ) which is given by the following lemma. 

Lemma 4.4 For all p £ N, and z £ Z, we obtain in the basis (e Q ) Q : 


= y J- 

V rv! 


M=p 


Proof. 


a*(z) p |fl) = a*(z 1 e 1 + ... + z K e K ) p |fi) = {z\a*{e{) + ... + z K a* (e K )) p \Ll) 


E ■ ■ ■ z K K **(e,r ■ ■ ■ a*(e K ) a * |n> 

|a|!=P 

p! 

a! 


= E ■ 

L ^ rv\ 

M —p 


And then 


a*{z) p ^ /p\ 


\/e p p\ 


|°>=£V& 


£ erv • 


|a|=p 


□ 

In the case where the initial state is a twin state, the following lemma is used to obtain an expansion in 
the basis (e a ). 


Lemma 4.5 Let <p , ij) be in Z and z the state 

a* (cf>) n a* (ip) r 


such that n + m = N. 
Then we obtain 


Vn!m! El ^ El 


Ve" +ra n!m! 

Vyl 


l«> 


| 7 |=iV |o:|=n,CK <7 


cd (7 — a)! 






Proof. 


E -E“ a *( e r 

ck! 

|q;|— n 

a*W m = E W^ a *( e ) fi 

\P\=m 

Wa*«o m = E 5ww a “‘(^= E 


| Q!|=7l, |/31 =771 

because of the CCR relations (1). 


|o|=7i,|/3|=ra 


n!m! 

cn!/3! 


a ,lfir,*(v\ a +P 


(j) a ip p a*{e) 


in) = V 

\fn\vi\ ^ old! /TTLdTmT 


| Q; | = 71 , |/31 —in 


a!/3! yETp 


in) 


= V n!m! E] 


| 7 |=AT,|a|=n,o :<7 


vP ^ a ,/,- T - a s*(e) 7 

V 7 ! 


cd (7 — a)! 


‘TV’ 


in) ■ 
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By setting 7 = a + /3 


= Vn\rn\ 'y^ 

ItI =n 


E 


|a|=n,Q :<7 


a!(7 — a)l 


c/) a ^- a ) 


«*(e ) 7 

\/e N j\ 


\n). 


□ 

Further the limit reduced density matrices ( 6 ) have to be computed numerically. In order to do this, the 
integration over Z of the Wigner measure is discretized. The problem is then reduced to the computation 
of the matrix elements of in the basis (e a )| Q | =p . 


Compute the matrix elements of |z® p )(z® p |, according to Lemma 4.4: 


*(e) c 


\Ze p a\ 


l«) = 



Then 


(|z® p )<z® p |)e a = Y, \ 

v m=p v p ' 


= E 

\p\=p 


-J== . 2 a Z% 


For the computation of the integral j z \z® p ) {z® p \dii(z), the Wigner measure is approximated by a convex 
combination of gauge invariant delta functions S z , where S z = f^ 5 e ie z d,0. 

For the Wigner measure associated with the Hermite states, n = S z 1 , and the discretization is trivial and 
exact, it is not needed to be approximated because of the gauge invariance. 

In the case of the twin states given by 41 n = a fe) la ^ 2 ' ) 2 |0), where ipi ,ip 2 € Z, ||^>i|| = ll^ll = 1, and 

•y e Tl 1 + ri 2n 1 !n 2 ! 

n\ = n 2 = y i the Wigner measure is fj,o = ^ 5^ d<j> according to [6] , with: 

VV = cos( 0 )^ o + sin, 
y/2 y/2 

4 >o = -yWh + ^2), V’f = i—Ol - ^2) • 

Numerically the interval [0, 27 t] is discretized and /io is approximated by <5f fc . 

The Wigner measure fi t propagated at the time t of the twin states is given by : 

2 /*27T ^ 

P ‘= 2 ni 


where is solution to the Hartree equation at the time t with initial condition 


Numerically it is now approximated by Z Y^k—i $z k (t) w here z k (t) solves the Hartree equation (7). 
Thus the matrix elements of J z \z® p ) (z® p \diit.(z) are given by the formula: 


And the approximation 

The matrix of 7 p (t) 
approximation. 


of 


P- 


1 m 

ly 
171 iti 

the scalar J z \z\ 2 p dfio(z) is 

— r 1 n— can then 

Jz \z\ 2 p dpo(z) 


z k (t) a z k (tf . 

given by the formula 
be computed at any 


time t numerically with a good 
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5 Error estimates 

5.1 Error estimate of the composition method 

The Baker-Campbell-Hausdorff formula (see [35]) allows to find the order of the composition method which 
is 4. Then the Taylor’s formula with 4* ,! order integral remainder and the Cauchy inequalities are used to 
estimate the error. 

The following proposition gives an estimate of the composition method. 


Proposition 5.1 Let R > 0, A and B be two anti-adjoint matrices such that (ai — 0 , 2 )||^4|| — ^§ a ||.B|| < R. 


Then 


\e A+B - v^bII < ( ( 0l - a 2 )\\A\\ - ^||B| 


3a 2 , 

T 


where 4 >a,b = e 2 B e aiA e 2 B e 2 B e a ' lA e 2 B e 2 B e aiA e 2 B is the composition method. 

Proof. 


S A + B _ <^ A B = e A+B — e 2 B e a ^ A e a i B e a i B e a 2 A e a i B e a i B e a i A e a j B 

— C A +B 0 a TBa 1 A 0 ai X a 2 B A o. 2 +^ B . a± B 


= \\e~ a - 1 ^ B e-^ A e-^ B e A+B - e a * A e^ B e a ' A e^ B \ 


Then for z £ C 

\\e< A+B ^* zA , zB \\ 

< e -^4 5a |/m(*)|||B|| e a 1 |/m(z)|||A|| e ^|/rn(z)|||B|| e |/m( Z! )|(||A|| + ||B||) 

+ e - a 2 \Irn(z)\\\A\\ e -^^\Im(z)\\\B\\ eai \Im(z)\\\A\\ e ^\Irn(z)\\\B\\ 

= e - ; ^|/m(z)|||B|| e a 1 |/m(z)|||A|| e iy-|/m(z)|||B|| ( ' e |/m(z)K||A|| + ||B||) +e -a 2 |/m(z)||| J 4||) 
= e |/m(z)|(-^||B||+o 1 ||A||)/ e |Zm(z)|(|| J 4||+||B||) + e -o 2 |7m(z)|||A|K 

< 2 e h™C)l(-¥l|S||+ai||A||) e -a 2 |/m(z)|(||A|| + ||B||) 


Let us consider the holomorphic function on C defined by: 

• 

Since the composition method is of 4 th order then for A £ K, the Taylor’s formula with integral remainder 
yields: 


Ia,b{ A) = f — 
j 0 

||/a,b(A)|| < ^ sup ||/i 5) s (t)|| • 
51 te[o,A] 
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By the Cauchy’s integral formula, we know for each t £ [0,1]: 


/£*(*) 


— / 

5! 


2*7T J\ z . 


<ir 

27r Jo 


< sup 

|/m( 2 :)|<l 


< 

sup 


Sa,b{z) 


dz 


302 | 


Hence for all Ar and Br such that (ai — a. 2 ) || >1/? || — ^H-BrII < R we obtain: 

||/^, BH (A)||<2A 5 e«,ifA<l. 

Let A and B be such that (ai — a 2 )||H|| — < R. 

By setting A R = - - — 3 ., IIBII and Br = - - 3a , IIRII , we obtain 

(a 1 -a 2 )\\A\\-^\\B\\ (ai-a 2 )|| A\\ -||B|| 


A = 


(a 1 -a 2 )||A||-^||H| 


R 


Ar , B = 


{ai ^a 2 ) \\A\\-^\\B\ 

R 


Bf 


and (on - a 2 )\\AR\\ - ^f-\\B R \\ = R. 
Then 


Ia,b{ 1 ) = fA R ,B R 

II/a,b( 1)|| < 2e 


' ( ni -o 2 )M||-^||Bi 

R 


HI ( a i-a 2 )\\A\\-^\\B\ 


R 


<2—(( ai -a 2 )\\A\\- 


3 a 2 


B I 


□ 


5.2 Error estimate of the approximated composition method 

The composition method is approximated by replacing e - l ^f dT (- A K) by its 4 th order Taylor expansion, with 
some normalization factor. 

Errors estimates for this modified composition method rely on the two following lemmas. 

Lemma 5.2 Let E be a normed vector space, J £ N*, (fj)j and (gj)j two maps sequences from E to E such 
that for all j £ {1, , J}: 

• fj is linear. 

• \\fj(u)\\ = IlSihdll = IMI for all u £ E. 

• \/u£E ||u|| < Q => \\fj(u) - gj{u)\\ < 6 

For uq £ E, ||uo|| < g set Uj = fj(uj- 1 ) and Vj = gj(vj- 1 ) with Vo = uq. 

Then we deduce \\uj — vj\\ < J6 . 


15 



Proof. Let us proceed by induction on J. 

For J = 0, 11Mo — Vo || = 0 < 0<5. 

Let us assume ||itj — uj|| < J5 with the hypotheses fullfilled for j € {1,_, J + 1}. 


uj +1 - v J+ 1 = f j+ 1 (uj) - g j+l (vj) 

= fj+i(uj) - fj+i(vj) + fj+i{vj) - gj+i(vj) . 

Since fj+\ is linear and unitary we obtain: 

||/ j+ i(uj) - f j+ 1 (vj)\\ = \\fj+i(uj - vj)\\ = ||wj - f j|| < JS . 

Moreover ||^j|| = \\gj o g J _ 1 o ... o 5l (w 0 )|| = ||u 0 || < g , 

then 

\\fj+i(vj) -g j+ 1 (vj)\\ < 6 , 

we obtain 

||«j+i - uj+i|| < (J + l)^ • 

□ 

Let TL{e A ) denote the A th order Taylor expansion of e A around 0. 

Let A be in B(0,cr), cr > 0. 

Lemma 5.3 Let u be a^vector in a normed vector space E and let A be an anti-adjoint operator on E. 
Define the application TL(e A ) on E which is non linear by: 

TL(e A )u = llT ^ M TL{e A )u if \\TL(e A )u\\ ? 0, 

it preserves the norm. 

Then ||(TL(e A )-TL(e A ))u|| < ||TL(e A ) - e A ||||u||. 

Proof. 


TL(e A )u - TL(e A )u = TL{e A )u - 
-(i- 


TL(e 


|| TL{e A )u 
U " ' TL{e A )i 


\\TL(e A )u\ 


\\TL(e A )u - TL{e A )u\\ = 


1 - 


||TL(e A )w 
= \\\TL(e A )u\\-\\u\\\ 
= \\\ TL (e A )u\\ - || e A u 
< || TL(e A )u — e A u\\ 

\\(TL(e A )-TL(e A ))u\\ < || TL(e A ) - e A ||H|. 


TL(e A )u 


□ 
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The 4 th order error of the Taylor expansion gives: 


— TL(e 


°° a 1 

Ev 


2=5 


<\\A\\ 5 r^^\\e tA \\dt<\\A\\ 5 ^. 


The following proposition gives an estimate of the approximation of the composition method. 


Proposition 5.4 Let A and B be two anti-adjoint matrices and J an integer such that 

^((a 1 -a 2 )\\A\\- 3 f\\B\\)<5 


and 


j>l_ {iai _ a 2 m \- 3 f\\B\\). 


Then 


II e^ A+B \ {*^ A ^ B ) J u\\ < - a 2 )||A|| - ^HBIl) 5 + |||^|| 5 ^ t^-\\u \| 

for all vector u, where 

q AB = e^TL(e aiA )e^e^TL(e a2A )e^e^TL(e aiA )e'^ L . 

Proof. 

Let u be a normed vector. 

First for i = 1,2,3 let us estimate the error: 

He - =— e aiA e~^~u - e~^~TL(e aiA )e~^~ u|| 

a i B 

by using the fact is an unitary operator and Lemma 5.3: 

||e^e° iA e^u - e^TL(e a ’ A )e^%|| = ||(e a ' A - TL{e aiA ))e^ u\\ 

< IMI | 5 < II^H 5 ^ 


60 


like in the previous proof. 

a * R A a i R a i B '—'— a a i B 

Now Lemma 5.2 can be applied with fi = e T e ai e T , gi = e V TL(e ai )eT and J = 3. 
Then 


ItE 'a,bu — H/a,bw|| < i l|A|| 5 . 


Secondly let us estimate the error: 


I e A+B u - 4 >a,bu\\ 


By using Proposition 5.1 with its hypotheses and the previous estimates: 

||e A+B u - ^a,bu\\ < \\e A+B - Ta,b||||m|| + ||(^a,b - $a,bH| 

-l?( (ai " a2)l|A|l " 3 f 2||B|l ) 5 + ? MI15 - 
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By applying that to and ^ B where At = j with J positive integer, we obtain: 


3 




< (I? ((ar-a 2 )P|| 



+ 4PII 5 


At 5 
£ 5 


Then by applying Lemma 5.2 with fj = and gj = $a t A a* b , we obtain: 


\e^ A+B) u (*^ A ^ B ) J u\\ < ( ^ (fa - a 2 )\\A\\ ^ ||B|| V + l\\A\A 


By knowing that for all positive integer r and positive real a, R 1 —> 
minfl >0 = (^r) T , the condition 


R T 


is minimal in R m i n 


^(V-a 2 )P||-^||B||) <5, 

that is 

J>^((a 1 -a 2 )||A||- 3 | 2 ||B||), 

implies 

\\e^ A+B) u &^ a> m b ) J u\\ < ^2 (|) 5 ((a, - o 2 )||A|| - ^||B||) 5 + t^- . 


- and 

a 


For Q £ C (\/ 2 Z), we know according to (15): 


Q 


Wick 

IV" 2 


Jn\(n + 2 -2)! 2+2 ~ 

- - -—t- e 2 5 „_ 2 + 2 (Q(g) 

(n — 2 )! 

= e 2 n{n - l)S , n (Q 18 ) Id® n ~ 2 ) , 


□ 


then 

\\QYv» z \\ < ^ N ( N ~ 1)11^11 IIQII < e 2 N 2 \\Q\\ = IIQH ■ 

When Q wick = V with 

V(ei V ej) — 2^1^^ ^ ; 

the norm ||V|| is bounded from above by ||V|| = \ max|Vjj| independently of the number N = ( 7 ) of 
particles. 


Moreover ||c?r( J 4)| yw z \\ < £AT||^4|| = ||A|| , therefore ||rfF(— Ax)|| < ||A^-|| = 2. 

Finally by applying the last proposition with A = — idr(— Ak) and B = — iV, an error estimate is obtained 
for the complete evolution: 




U ( V I'_ At irf r(-A*r),- &± iv ) J u\\ < 



a 2 ) II Ajc|| 





Pratically, the time step is chosen according to N and t so that the above error is negligeable. 
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6 Numerical simulations 

For all the numerical simulations the final time is chosen to be t max = 1, the number of time steps for the 
4 th order Runge-Kutta method applied to solve the mean field equation is 100. 

The loop of the number of particles is performed numerically from N m i n = 2 to N max = 20 particles, and 
only for an even number of particles. 

In the Fortran program the computations were performed by parallelizing the loop in the computation of 
the product sparse matrix-vector dr(— Ak)u with Openmp on 8 threads. 

Results and orders of convergence for and y| 2 ' 

For each type of states, the following graphics show for the reduced density matrices and for K = 10 
sites: 

1) The logarithm of the error in trace norm log(max tg [ 0 

of the number of particles N in the cases p = 1 and 2. 

A straight line is obtained whose the slope is the order of the error in 1/N. 

These numerical experiments also valid the idea that for rather smooth but non trivial iV-body bosonic 
system, the mean field asymptotics start to be relevant at N = 4. The numerical plot agree perfectly 
with the theoretical results of [ 1 ]. 

2) In the case p = 1 the density of particles on each site k £ {1,..., A'} given by 7 ^ (t) for N = 20 
particles and for the mean field limit at the same times t = 0 and t = \. 


7 ffit) — 7 ; x’{t) ) according to the logarithm 


„(p) t 


3) The correlations in terms of the 1 and 2 particles reduced density matrices, for N = 20 particles and 
the mean field at the time t = 1. Depending on the case, this plot shows with which accuracy the 
mean field also catches some quantum correlations. 

6.1 Hermite states 

For the Hermite states z® N the vector £ is given by 2 = ^((1 + i)ei + ie 3 ). 


Log{ max||7 j\*’(“ 7^, J (<)|| ) according to ,N € [2,20] , K = 10 , p = 1 



(a) Log-log plot for Hermite states. p=l, K=10. 
Numerical slope -0,9804 


Time evolved densities of particles on each sites k, given N and mean field 



k 


(b) Compared densities of particles at times t=0 and t=l 
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Log{ max Ikv ' ( t) — 7 i^(0|| ) according to Loq{N) . A r € [2.20] . K — 10 . p = 2 
te[o.i]ll 111 


Erroi=0.36 - 


- LOG ERROR IN TRACE NORM 

S LO P E :-0.9693643003528060742013 



Log(N) 

(c) Log-log plot for Hermite states. p=2, K=10. (d) Mean field(white) and 20-body quantum(blue) correlations 

Numerical slope -0,9693 at time t = 1 


6.2 Twin states 

For the twin states 'Fjv = ° la ^ 2 ' > ~ |0), ipi = -L(ei + ies) and tp 2 = ■ 

E n i + n 2 m !n 2 ! v2 


Loq{ max 

* \e[o.i] 


| 7 _^(t) — 7^(t)|| ) according to Log(N) ,N € [2,20] , 


K = 10, p = 1 


Erroi=0.135 - 


- LOG ERROR IN TRACE NORM 

SLOPE: -0.9855695390119798560136 


0.082- 


0.049- 


0.030 - 


0.018- 


0.011- 



ErroF 0.006 H- 1 -■- 1 -1- 1 - 1 -»-1- 1 -'- 1 -1- 1 - 1 -’-1- 1 -1- 1 -1 

0.5 1.0 1.5 2.0 2.5 3.0 

Log(N) 


Tiine evolved densities of particles on each sites k. given N and mean field 



k 


(e) Log-log plot for twin states. p=l, K=10. 
Numerical slope -0,9855 


(f) Compared densities of particles at times t=0 and t=l 
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Log( max || 7 y J (£) — )|| ) according to Log{N) ,N € [2,20] , K — 10 , p = 2 

- LOG ERROR IN TRACE NORM 



Log(N) 

(g) Log-log plot for twin states. p=2, K=10. 
Numerical slope -1,0566 


((Awj - (tSoM7»)jj)(1) ((7 Who - (AfrtthjHi) 



(h) Mean field(white) and 20-body quantum(blue) correla¬ 
tions at time t = 1 


6.3 Wq states 

For the Wq states 'F n 


a* (ifti) 711 a* (^ 2 ) 712 

yj e n i+ n 2m!n2! 


|fi), V’i 


= 75 ( ei + ie3 ^ and = e2 ‘ 


In this case the state is given by Sn(iI> f Ul <S> with n\ + ri 2 = N, n\ = N — q and n 2 

the mean field. 

The associated Wigner measure is . 

In these simulations q = 2. 


q fixed for 


Lo>A wax IhUV) — ill 1 (Oil) a«»ixlii« to Lag(S), N € [2.20], A'= I0.;» = 1 




(i) Log-log plot for Wq states. p=l, K=10. 
Numerical slope -0,9843 


(j) Compared densities of particles at times t=0 and t=l 
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E no 1=2.22 -, 


W wax |h2'(0- ^'(1)11) a«oixlii« to/.osr(A0,JV€ [2,20], A'= 10,j> = 2 



(k) Log-log plot for Wq states. p=2, K=10. 
Numerical slope -0,9449 


((Awj - bff MtS’JjjXI) ((l£ ] )ij,ij - (7^ ) )M(7il ) U(l) 



(1) Mean field (white) and 20-body quantum(blue) correlations 
at time t = 1 


6.4 Other states 


A case when the order of convergence is equal to l/2.(see[l]) 
In this case g e = \(/)% N )with <^v = -Lei + \J 1 - ^e 2 . 
The associated Wigner measure is 5^. 


Log{ max ||7jy ; (t) — 7 ^(<)|| ) according to Log(N) ,N € [2.20] , K = 10 , p = 1 


- LOO ERROR IN TRACE NORM 

SLOPE: -0.4711851936408901764253 



Time evolved densities of particles on each sites k. given N and mean field 



k 


(n) Compared densities of particles at times t=0 and t=l 
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Loq{ max 
<e[<u] 


Erroi=1.822 


1.401- 


1 . 221 ' 


0.818 


0.670 


Error=0.548 


JTjv^) ” 7 ^ , (£)|| ) according to Log(N) ,N € [2,20] , K = 10 , p = 2 


- LOO ERROR IN TRACE NORM 

SLOPE: -0.4624099403778836467150 



16 Log(N) 


(o) Log-log plot. p=2, K=10. Numerical slope -0,4524 (p) Mean field(white) and 20-body quantum(blue) correla¬ 

tions at time t = 1 
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Appendix 

Class of symbols [3, 4, 5, 6]. For any p, q £ N, define & P ,q to be the space of homogeneous complex¬ 
valued polynomials on iF such that b £ VA if and only if there exists a (unique) bounded operator 
b £ if(V p ^,V 9 JF) such that for all 2 £ S’: 

b(z) = (z® q ,bz ® p ). (14) 


uWick _ -i y/n\(n + q p)\ ^£±2 c ®{n—p)\ ocl 

0 \V n 3f — l[p,+oo)(^j ^ £ ^n—p+q (^0 ® 1 J , (lb) 

where 6 denotes the operator associated with the symbol b according to (14). 

The composition method based on the Strang splitting with the coefficients (12) is of 4 th order (see [35]). 


Dimension of \J N C K : ^ i ^ ^ or ^ an d N in [1,20]: 



K = 1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

N = 1 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

2 

1 

3 

6 

10 

15 

21 

28 

36 

45 

55 

3 

1 

4 

10 

20 

35 

56 

84 

120 

165 

220 

4 

1 

5 

15 

35 

70 

126 

210 

330 

495 

715 

5 

1 

6 

21 

56 

126 

252 

462 

792 

1287 

2002 

6 

1 

7 

28 

84 

210 

462 

924 

1716 

3003 

5005 

7 

1 

8 

36 

120 

330 

792 

1716 

3432 

6435 

11440 

8 

1 

9 

45 

165 

495 

1287 

3003 

6435 

12870 

24310 

9 

1 

10 

55 

220 

715 

2002 

5005 

11440 

24310 

48620 

10 

1 

11 

66 

286 

1001 

3003 

8008 

19448 

43758 

92378 

11 

1 

12 

78 

364 

1365 

4368 

12376 

31824 

75582 

167960 

12 

1 

13 

91 

455 

1820 

6188 

18564 

50388 

125970 

293930 

13 

1 

14 

105 

560 

2380 

8568 

27132 

77520 

203490 

497420 

14 

1 

15 

120 

680 

3060 

11628 

38760 

116280 

319770 

817190 

15 

1 

16 

136 

816 

3876 

15504 

54264 

170544 

490314 

1307504 

16 

1 

17 

153 

969 

4845 

20349 

74613 

245157 

735471 

2042975 

17 

1 

18 

171 

1140 

5985 

26334 

100947 

346104 

1081575 

3124550 

18 

1 

19 

190 

1330 

7315 

33649 

134596 

480700 

1562275 

4686825 

19 

1 

20 

210 

1540 

8855 

42504 

177100 

657800 

2220075 

6906900 

20 

1 

21 

231 

1771 

10626 

53130 

230230 

888030 

3108105 

10015005 
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